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ABSTRACT 

Recent  experiments  on  the  TMX-U  device  at  LLNL  have  indicated  the 
possibility  of  a  drift  wave  being  driven  unstable  by  the  injection  of  neutral 
beams  in  the  thermal  barrier  region.  A  review  of  linear  theory  is  presented  in 
the  local  approximation.  In  addition,  particle  simulations  are  used  to  understand 
the  nonlinear  characteristics  of  this  instability.  The  particle  simulations  are  per¬ 
formed  using  Id  -  3v  and  2d  -  3v  electrostatic  particle  simulation  codes.  The 
instability  is  shown  to  saturate  by  beam  trapping,  nonlinear  electron  effects,  and 
nonlinear  motion  parallel  to  the  density  gradient.  Resonant  ExB  motion  is  a 
significant  process  for  some  of  the  beam  particles.  This  process  is  closely  asso- 
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dated  with  the  trapping  of  beam  ions.  Harmonic  generation  is  also  observed,  as 
appears  in  some  of  the  experimental  data.  Qualitative  nonlinear  theory  is 
presented  in  support  of  the  particle  simulations. 
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1.  Introduction 

The  theoretical  study  of  collisionless  drift  waves  was  begun  a  number  of  years  ago1.  In 
particular,  stability  boundaries  were  given  for  the  electrostatic  electron  drift  wave  in  that  origi¬ 
nal  paper.  The  instability  mechanism  given  in  that  paper  was  inverse  Landau  damping  on  the 
Maxwellian  electron  component.  The  Maxwellian  ion  component  served  to  Landau  damp  the 
waves.  The  overall  stability  was  then  governed  by  the  difference  in  the  effects  from  the 
resonant  electrons  and  resonant  ions.  Particle  simulations  have  also  been  performed  for  this 
instability  as  in  Ref.  2  and  Ref.  3. 

In  this  paper  the  analysis  is  extended  to  include  the  effects  of  injected  ion  beams  as  might 
occur  in  the  thermal  barrier  cells  of  a  tandem  mirror  device  .  For  this  case  there  is  extra  free 
energy  associated  with  the  ions  which  can  cause  instability  .  This  ion-beam-driven  drift 
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instability  is  a  likely  candidate  for  explaining  some  of  the  early  experimental  results  on  TMX-U 
at  LLNL.  It  is  in  this  context  that  the  instability  is  of  practical  interest.  This  paper  concentrates 
upon  the  instability  for  parameters  similar  to  that  of  present  day  thermal  barrier  cells. 

Section  2  reviews  the  linear  theory  using  the  electrostatic  local  approximation.  Of  partic¬ 
ular  interest  is  the  dependence  of  growth  rates  and  wavelengths  on  the  relative  density  of  the 
beam  to  the  background  plasma,  the  thermal  spread  of  the  beam  along  the  magnetic  field,  and 
the  zero  order  beam  velocity  perpendicular  to  the  magnetic  field. 

Section  3  presents  the  simulation  models  and  linear  simulation  results.  Section  4  includes 
the  qualitative  nonlinear  theory  and  the  nonlinear  results  from  the  simulations.  The  characteris¬ 
tic  nonlinear  behavior  includes  nonlinear  electron  effects,  mode  coupling,  ion  beam  trapping, 
and  enhanced  cross  field  motion  of  trapped  beam  ions.  Section  5  includes  a  short  discussion  of 
the  effects  of  axial  inhomogeneity  on  the  linear  and  nonlinear  characteristics  of  the  instability. 
A  summary  and  conclusions  are  given  in  Section  6.  Appendix  A.  presents  a  discussion  of  the 
field  solves  used  in  the  particle  simulation  codes. 

2  Linear  Theory 

In  this  section  linear  theory  is  examined  in  slab  geometry  using  the  local  approximation. 
The  effects  of  the  various  parameters  of  the  problem  are  noted  and  the  most  unstable 
configurations  are  identified.  Stabilizing  effects  are  also  presented. 

A.  Theoretical  Model 

The  frequency  regime  of  interest  for  this  problem  is  o>  «  Clci.  In  addition,  it  is  assumed 
that  kj)e  «  1  and  c o/k.  «  vle.  Only  the  electrostatic  response  is  kept  which  is  valid  if  the 
ion  beam  velocity  Ub  (  along  z  )  is  much  lower  than  the  Alfven  velocity.  Variation  of  only  the 
plasma  density  is  considered  and  the  magnetic  field  is  assumed  to  be  uniform  in  space.  Under 
these  conditions  the  dispersion  relation  may  be  written4 

Z)(k,o,x)  •=  =  0  (1) 
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<»p j  is  the  unperturbed  plasma  frequency  evaluated  at  a  location  x0,  ex  =  (1/2)  vf ,  f,  (A’,€i,v.) 
is  the  unperturbed  guiding  center  distribution  function  with  X  =  x+vylClc  ,  the  argument  of 
the  Bessel  function  is  kj>s,  and 
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The  plasma  inhomogeneities  are  in  the  x  direction  which  is  perpendicular  to  the  uniform  mag¬ 
netic  field  in  the  z  direction.  For  consistency  with  the  local  approximation  it  is  required  that 
(vJfl^d/idX)  «  1  and  that  ky  »  d/(dX). 

For  this  parameter  regime  the  electron  response  reduces  to  Debye  shielding  The  ion 
beam  and  background  ion  susceptibilities  are  more  complicated.  The  distribution  function  for 
the  beam  or  background  ions  may  be  taken  as 
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where  a  is  a  constant,  v„  =  y/lT,/m,)  being  the  ion  parallel  thermal  speed,  v0  =  0  for  the  back¬ 
ground  ions  and  v0  =  ±  Ub  for  the  counter  streaming  beam  ions.  This  model  allows  finite  kj), 
effects  and  thermal  effects.  The  model  also  makes  it  possible  to  consider  the  case  where  the 
zero  order  vx  and  v.  are  comparable,  which  is  a  case  of  considerable  practical  interest. 

For  the  moment,  however,  consideration  will  be  given  to  the  case  where  the  ion  beams 
and  the  background  ions  are  completely  cold  and  the  beam  ions  have  zero  gyroradius.  In  this 
cold  limit  the  physics  is  especially  clear  and  the  dispersion  relation  becomes 
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where  8  is  the  ratio  of  the  beam  density  to  the  total  ion  density,  and 
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represents  the  inverse  scale  length  of  the  background  ion  population  evaluated  at  the  location 
x  *  x0  where  the  local  approximation  is  being  employed.  (  Kb  *■  1/Lb  represents  the  inverse 
scale  length  in  the  ion  beam  density)  The  terms  proportional  to  8  are  from  the  ion  beams.  The 
terms  with  the  frequency  dependence  of  \/o)2  represent  the  parallel  inertia]  response,  the  terms 
with  the  frequency  dependence  of  \/<d  represent  the  ExB  response.  The  other  ion  derived  term 
in  the  equation  is  from  the  polarization  drift,  which  in  this  approximation  of  zero  gyroradius  is 
the  same  for  both  species  of  ions.  Finally,  only  Debye  shielding  is  retained  for  the  electrons. 

When  8  is  an  infinitesimal  quantity  and  the  beam  spatial  gradients  may  be  neglected  the 
solutions  to  Eq  (5)  are  the  drift  wave  branch 
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Here  cs  is  the  sound  speed  and  ps  =  cs/Clci.  The  solutions  given  in  Eq.  (7)  and  .Eq.  (8) 
assume  that  the  system  is  quasineutral  The  solutions  given  in  Eq.  (7)  and  Eq.  (8)  are  shown  in 
Fig.  1.  Notice  that  when  Ub  >  cs  and  Ub  <  cs/(k.2Ln)  the  beam  modes  and  the  fast  branch 
of  the  drift  wave  intersect,  indicating  that  the  the  two  modes  are  strongly  coupled  and  instabil¬ 
ity  may  result.  In  Fig.  2.  a  numerical  solution  to  the  dispersion  relation  for  a  beam  with  finite 
density.  Only  the  real  portion  of  the  frequency  is  plotted. 


The  approximate  scaling  for  the  parameters  of  this  system  can  be  obtained  from  Fig  1. 
Specifically,  a >r  =  Ubk:,  kj>s  =  1,  and  k.L„ ^  cs/(2Ub).  The  quasineutral  approximation  is 
good  in  the  limit  ftf/  «  cdp,  ,  which  corresponds  to  the  case  of  interest  for  thermal  barrier 
cells.  Note  that  the  value  for  k^L,,  remains  unspecified.  For  consistency  with  the  local  approxi¬ 
mation  it  is  required  that  kJLn  »  1,  which  is  valid  in  the  limit  a >  «  flf/. 


Note  that  the  usual  field-aligned  electrostatic  ion-ion  two-stream  instability  is  unstable 
under  the  opposite  conditions  in  the  cold  limit;  namely  Ub  <  cs.  Since  the  drift  wave  is  merely 
an  extension  of  the  sound  wave  for  a  situation  where  there  is  a  density  gradient,  there  is  a  close 
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connection  between  these  instabilities.  In  some  limits  there  will  not  be  a  unique  way  to  name 
the  instability  since  it  may  exist  both  with  and  without  the  density  gradient. 

B.  Variation  of  Parameters 

Variation  of  the  characteristic  parameters  is  expected  to  change  the  situation  considerably 
from  the  cold  limit.  Some  parameter  variations  are  considered  in  this  section.  The  most  impor¬ 
tant  of  these  parameters  are  the  thermal  spread  of  the  ion  beam  parallel  to  the  magnetic  field  , 
finite  gyroradius  effects,  the  beam  fraction,  the  electron  temperature,  and  the  density  gradients 
(  which  need  not  be  the  same  for  the  different  species  ).  Attempting  to  vary  all  of  the  parame¬ 
ters  is  intractable  and  so  only  variation  of  one  beam  parameter  at  a  time  is  considered.  In 
addition,  the  background  component  is  always  assumed  to  be  cold. 

The  plots  of  real  and  imaginary  frequency  will  be  scaled  to  the  quantity  a^ax  *  . 5cs/Ln , 
which  frequency  represents  the  maximum  drift  frequency  for  cold  ions  and  vanishing  kz .  It 
should  be  remembered  however,  that  the  analysis  presented  in  this  paper  is  valid  only  for 

cu  «  nc/. 


The  general  susceptibility  for  the  ion  beams  including  all  of  the  parameters  of  the  model 
is  given  by 
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where  Z  is  the  usual  plasma  dispersion  function  ,  ATf  =  f-,  £  «*  ^xVi0 ,  and 
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aj,  =  (£iV„2)/(  ft  „,£,,,).  The  set  of  parameters  that  are  being  used  as  a  reference  case  are  given 
in  Fig.  2  and  will  be  fixed  unless  otherwise  noted. 


Starting  with  a  given  set  of  parameters,  it  is  desired  to  vary  a  given  parameter  and  to  fol¬ 
low  the  most  unstable  root.  The  equations  that  must  be  solved  to  find  the  most  unstable  root 
are 
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where  k  is  taken  to  be  a  real  quantity  and  the  equations  are  solved  for  complex  w.  The  two 
equations  given  in  Eq.  (10)  and  Eq(ll)  constitute  four  nonlinear  equations  for  the  four  unk¬ 
nowns  Aj_ ,  k.  ,  the  real  part  of  the  frequency  a>r  and  the  growth  rate  y.  Efficient  mathematical 
routines  exist  for  solving  such  systems  of  equations5  These  routines  are  used  in  a  computer 
program  written  by  one  of  us  (  Y.-J.  Chen  )  and  were  very  valuable  in  determining  the 
behavior  of  the  linear  dispersion  relation.  Figures  (3)- (7)  were  generated  by  this  computer  pro¬ 
gram. 

One  of  the  most  striking  effects  is  shown  in  Fig.  3,  which  demonstrates  the  effect  of  hav¬ 
ing  the  ion  beam  gradient  not  equal  to  the  gradient  of  the  background  plasma.  In  particular,  the 
growth  rate  is  highest  when  the  spatial  gradients  for  the  two  species  are  in  opposite  directions. 
This  is  not  a  small  effect.  This  effect  is  also  possible  to  realize  in  experiments  since  the  ion 
beam  is  created  by  external  sources.  For  example  when  the  optical  depth  of  the  neutral  beam  is 
short  compared  to  the  width  of  the  system,  the  density  gradients  for  the  two  ion  species  will  be 
in  opposite  direction.  In  all  of  the  following  plots  the  ion  beam  gradient  is  set  equal  to  zero  so 
as  to  give  an  estimate  of  the  more  unstable  situations  which  may  occur.  The  qualitative  effects 
of  changing  the  other  parameters  are  the  same  regardless  of  the  value  of  the  beam  density  gra¬ 
dient. 

The  effect  of  finite  kj>,  for  the  beam  ions  is  most  clearly  seen  from  the  susceptibility 
given  in  Eq.  (9).  There  it  is  seen  that  the  destabilizing  parallel  drift  terms  are  multiplied  by 
y02  (kj>,)  and  so  as  kj),  becomes  on  the  order  of  unity  the  effective  beam  density  decreases. 
Since  could  be  at  least  comparable  to  as  Ub,  the  interesting  values  of  kj>,  are  could  also 
be  of  order  unity.  A  more  complete  description  of  finite  gyroradius  effects  is  given  in  Fig  4  for 
a  typical  example.  As  the  gyroradius  of  the  beam  ions  is  increased  the  growth  rate  decreases, 
the  value  of  kj>s  decreases,  and  the  real  frequency  also  increases  All  of  these  effects  are  con 
sistent. 

Variation  of  the  beam  parallel  energy  spread  also  reduces  the  growth  rate  of  the  most 
unstable  mode.  This  is  seen  in  Fig.  5.  where  the  most  unstable  mode  for  a  given  case  is 
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followed  as  the  thermal  spread  is  increased.  The  presence  of  resonant  ions  also  forces  the  phase 
velocity  of  the  wave  along  the  field  lines  to  decrease.  The  growth  rate  does  not  become  nega¬ 
tive,  however.  It  is  worth  mentioning  that  the  thermal  spread  is  a  slight  destabilizing  factor  for 
one  root  which  is  otherwise  stable. 


In  order  to  understand  the  electron  effects  more  completely  it  is  necessary  to  include  the 
imaginary  portion  of  the  electron  susceptibility.  This  is  done  by  using 
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where  Z  is  the  usual  plasma  dispersion  function  ,  ae  =  Vl2T./me),and  <o'= - k„  with 
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Ke  the  inverse  scale  length  of  the  electron  population.  Consider  the  addition  of  a  small  popula¬ 
tion  of  cool  electrons  For  these  electrons  to  >  c o‘00,  and  the  most  unstable  mode  is  stabilized. 
Again  as  in  the  case  of  the  ion  thermal  spread  a  mode  which  is  stable  with  no  cool  electrons 
may  made  slightly  unstable  by  the  addition  of  the  cool  electrons. 


Now  consider  the  electrons  consisting  of  a  single  hot  species.  Then  a>  <  w'  and  the  most 
unstable  mode  is  slightly  destabilized.  This  effect  is  one  effect  which  may  cause  the  destabiliza¬ 
tion  of  drift  waves  with  no  free  energy  in  the  ion  velocity  space.  The  effect  is,  however,  very 
small  when  (o/(k.ae)  «  1  which  corresponds  to  the  cases  of  interest  here.  It  should  be  noted 
that  this  discussion  of  the  electron  effects  has  assumed  that  the  electron  density  gradient  is  in 
the  same  direction  as  that  of  the  background  ions.  This  assumption  need  not  be  true  for  the 
energetic  electron  component  since  their  spatial  distribution  is  greatly  influenced  by  external 
factors  (  such  as  external  heating  ).  Therefore  the  electrons  may  serve  to  increase  or  decrease 
the  growth  rate  of  the  most  unstable  mode,  but  in  general  cause  only  a  small  correction. 

Figure  6  shows  the  effect  of  variation  of  the  relative  beam  to  background  density.  Notice 
that  the  maximum  growth  rate  occurs  for  a  value  of  8  less  than  one  half.  Also  notice  the  other 
variables  presented  in  Fig.  6.  First  of  all  ,  as  the  beam  fraction  approaches  zero  the  real  fre¬ 
quency  approaches  Ubk.  as  is  expected.  In  this  weak  beam  limit  the  growth  rate  y  is  propor¬ 
tional  to  the  square  root  of  the  beam  density  and  the  growth  rate  does  not  become  zero  (  in  the 
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cold  limit  )  until  the  beam  density  itself  becomes  zero.  In  this  weak  beam  limit  the  quantity 
kzLn  has  its  maximum  value. 

As  the  beam  fraction  increases  the  real  frequency  becomes  lower  for  two  reasons.  First 
the  beam  modes  begin  to  be  split  from  the  bounce  frequency  as  in  Fig.  2.  Secondly,  the  max¬ 
imum  drift  frequency  propagating  on  the  background  plasma  becomes  smaller  as  the  back¬ 
ground  density  becomes  smaller  Thus  there  would  be  no  intersection  of  the  beam  modes  with 
the  drift  wave  branch  unless  the  frequency  of  the  beam  mode  was  lowered.  This  is  clearly  seen 
in  Fig.  6.  where  the  real  frequency,  the  growth  rate  and  the  value  for  kzLn  are  all  decreasing. 

As  the  background  beam  density  goes  to  zero,  the  interaction  between  the  beam  modes 
themselves  needs  to  be  considered  since  the  beam  modes  may  have  some  of  the  properties  of 
the  drift  waves.  When  the  background  ion  density  does  go  to  zero  and  the  ion  beams  have 
their  density  gradients  in  the  same  direction  (  as  in  our  case  ),  no  instability  is  present  in  any 
parameter  regime. 

Finally,  the  quantity  Ub/cs  is  allowed  to  vary  and  the  results  are  shown  in  Fig.  7.  The  ins¬ 
tability  becomes  more  and  more  severe  as  the  streaming  velocity  approaches  the  sound  speed. 
As  this  happens  the  real  frequency  and  k.  both  increase  very  rapidly  and  the  value  for  kjys 
goes  to  zero.  In  this  limit  the  instability  reduces  to  the  field-aligned  ion-ion  two-stream  instabil¬ 
ity6'"  which  has  frequencies  and  growth  rates  on  the  order  of  o>p,  and  the  original  approximation 
of  keeping  only  the  low  frequency  response  is  not  valid.  For  those  cases  where  the  beam  velo¬ 
city  is  very  close  to  the  sound  velocity  the  temperature  of  the  beam  and  of  the  background 
become  very  important.  Severe  instabilities  for  cold  beams  may  prove  to  be  stable  for  more 
realistic  distribution  functions7. 

3.  Simulation  Model  and  Linear  Results 

This  section  describes  the  simulation  techniques  and  the  linear  behavior  of  the  instability 
as  recovered  from  our  particle  simulations. 


A.  Simulation  Model 
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The  simulation  geometry  is  presented  in  Fig.  8.  The  uniformity  of  the  magnetic  field  is  a 
consequence  of  the  low  beta  approximation.  As  in  Sec.  2,  the  electrostatic  approximation  is 
used  and  is  valid  if  the  beam  velocity  Ub  is  much  less  than  the  Alfven  speed  VA .  The  simula¬ 
tion  codes  use  nonlinear  Boltzmann  electrons  and  particle  ions  to  solve  Poissons  equation  (  or 
quasineutrality  )  in  the  y-z  plane  with  a  constant  magnetic  field  in  the  z  direction.  The 
quasineutral  approximation  is  valid  in  the  limit  w2,  »  ft2  ,  which  is  a  realistic  limit.  The 
time  step  constraint  is  determined  by  stability  of  ion  cyclotron  waves.  However,  this  does  not 
mean  that  «  1  if  only  the  particle  guiding  center  motion  is  desired.  The  stability  con¬ 

straint  due  to  the  ion  cyclotron  waves  can  be  removed  by  schemes  using  a  biased,  first  order 
accurate  time  integration  scheme8,  although  this  was  not  done  for  this  paper.  For  the  1-d  3-v 
code  only  one  direction  is  kept  in  the  y—z  plane  and  periodic  boundary  conditions  are  used  for 
both  the  particles  and  the  fields.  The  boundary  conditions  applied  to  the  fields  for  the  2-d  3-v 
code  are  periodic  in  the  y  direction  and  E„  =  0  at  z  =  0  and  z  ~L. .  The  particles  are  periodic 
in  the  y  direction  and  specularly  reflected  in  the  z  direction.  Therefore  this  2-d  3-v  simulation 
model  is  a  type  of  square  -  well  model  for  the  thermal  barrier  region.  Only  those  Fourier 
modes  which  are  important  for  the  instability  are  kept  self-consistently  in  the  codes. 

The  density  variation  is  in  the  x  direction  and  is  implemented  using  the  "  ghost  particle  " 
technique  from  Cohen9.  The  ghost  particle  technique  assumes  a  fixed  density  profile  for  the 
various  species  and  has  the  drawback  that  the  density  profile  cannot  be  self-consistently 
updated.  Thus,  final  saturation  levels  from  the  simulations  can  be  expected  to  retain  only  quali¬ 
tative  information.  This  ghost  particle  technique  is  useful  when  velocity-space  relaxation  may 
saturate  the  instability,  as  for  the  instability  being  studied  in  this  paper.  The  exact  form  of  the 
density  variation  may  be  specified  arbitrarily.  For  our  simulations  we  have  found  it  convenient 
to  use  rt(x)  —  n0(  1  +  rx/(l  +  lrx]))  where  rx  -  x/Ln  is  the  displacement  relative  to  the  local 
density  scale  length.  In  order  to  restrict  the  parameter  space,  the  beam  scale  length  was  taken 
to  be  infinite,  which  represents  the  more  unstable  possibilities  as  seen  in  Fig.  3.  Since  the  ions 
must  have  a  x  coordinate  associated  with  them,  the  2-d  3-v  code  has  six  phase  space  variables 
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and  the  1-d  3-v  code  has  five  phase  space  variables. 

B.  Linear  Simulation  Results 

In  this  section  representative  linear  code  behavior  is  presented.  Excellent  agreement  with 
theory  is  obtained  for  the  simplified  distribution  functions  used  in  the  simulations. 

Initially  we  consider  the  case  of  parallel  injection  only;  that  is  vx  «=  0.  A  series  of  test  runs 
was  conducted  in  order  to  compare  the  theoretical  growth  rates  and  those  given  by  the  2-d  3-v 
simulation  code.  The  simulations  parameters  were  8  =  0.2,  Ub/cs  =  3  0,  fic,Ar  =  0.5  -  1.0, 
‘“max/flf,  *  .03  and  only  one  Fourier  mode  was  kept  The  normalized  parallel  drift  velocity, 
vJUb,  was  varied  and  the  system  lengths  were  changed  so  that  the  most  unstable  mode  was 
retained.  There  were  Ny  =  65  grid  cells  in  the  y  direction  and  N.  *=  32  grid  cells  in  the  z  direc¬ 
tion.  Approximately  25,000  particles  were  used.  For  these  parameters  the  kJL„  ^  10  and  the 
local  approximation  can  be  expected  to  retain  some  validity,  at  least  in  the  linear  regime.  The 
parallel  velocities  were  chosen  by  bit-reversing  technique  as  in  Birdsall  and  Langdon10  This 
allowed  much  cleaner  results  than  from  a  random  initialization. 

An  example  is  given  in  Fig.  9  of  the  square  of  the  mode  amplitude  as  a  function  of  time. 
Note  that  the  line  is  not  perfect.  This  is  to  be  expected  since  for  this  simulation  resonant  parti¬ 
cles  play  a  large  role  in  the  determination  of  the  growth  rate.  Also,  the  simulation  could  be 
started  without  initial  excitation  in  which  case  the  total  growth  in  field  amplitude  would  be 
about  double  what  it  is  here.  The  characteristics  of  the  results  are  the  same  when  this  is  done. 
Finally,  the  five  test  cases  are  given  in  Fig.  10  along  with  the  calculated  growth  rates  from  Sec. 
2.  The  boundary  conditions  used  by  the  2-d  3-v  code  do  not  change  the  growth  rate  from 
infinite  medium  theory.  Agreement  is  very  good.  The  error  is  on  the  order  of  the  error  made 
in  the  theoretical  treatment;  that  is  the  error  is  on  the  order  of  oi/Cla. 

The  injection  angle  of  the  neutral  beams  relative  to  the  magnetic  field  in  the  thermal  bar¬ 
rier  cell  is  designed  to  be  approximately  forty  five  degrees".  Therefore  the  zero  order  vx  should 
be  approximately  equal  to  the  zero  order  streaming  motion  along  the  magnetic  field.  This  can 
be  and  has  been  added  to  the  simulation  model.  However  the  additional  noise  associated  with 
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the  ion  cyclotron  waves  requires  a  very  large  number  of  particles  to  suppress.  Thus  the  particle 
simulations  become  prohibitively  expensive.  One  could  try  to  use  a  guiding  center  mover12  or  a 
gyrokinetic  approach13  in  order  to  eliminate  ion  cyclotron  modes  and  still  retain  finite  gyrora- 
dius  effects.  However,  because  the  instability  characteristics  are  insensitive  to  finite  gyroradius, 
finite  gyroradius  effects  will  be  ignored  for  numerical  efficacy.  Specifically,  the  linear  characteris¬ 
tics  are  insensitive  to  finite  kj>,  corresponding  to  the  most  unstable  mode.  Similarly,  quasil- 
inear  effects  are  only  modified  by  a  multiplicative  factor  Joikj},)  from  the  case  with  zero 
gyroradius4. 

4.  Nonlinear  Behavior 

This  section  describes  the  nonlinear  phenomena  as  observed  from  the  simulations  and 
presents  a  qualitative  discussion  of  their  importance.  The  nonlinear  phenomena  include  (1) 
nonlinear  electron  effects  ,  (2)  nonlinear  motion  in  the  density  gradient  direction,  (3)  mode 
coupling,  and  (4)  trapping  of  the  beam  ions.  The  relative  importance  of  these  effects  depends 
upon  the  (1)  the  linear  growth  rate  ,  (2)  the  parallel  thermal  spread  in  the  ion  beams  and  (3) 
the  relative  density  of  the  ion  beams  to  the  background  plasma. 

A.  Motion  in  a  Single  Wave 

The  unstable  modes  that  are  present  have  their  phase  velocities  lower  in  absolute  values 
than  the  beam  velocity  as  seen  in  Fig.  (2).  The  qualitative  condition  for  trapping  depends  upon 
whether  the  ion  beam  is  warm  or  cold.  Here  the  beam  would  be  labeled  as  warm  if  many  beam 
particles  were  resonant  with  the  parallel  phase  velocity  of  the  unstable  wave  If  the  beam  is 
thermal  in  that  sense,  then  saturation  could  be  due  to  the  nonlinearities  associated  with 
resonant  particles.  If  this  happens,  then  saturation  occurs  within  about  an  e-folding  time  of 
when  the  growth  rate  is  equal  to  the  bounce  frequency  of  the  resonant  particles  in  the  wave. 
Estimation  of  the  perturbed  potential  at  saturation  by  this  method  yields  the  scaling  law 


(13) 


where  Ub  -  ac,  and  oib  =  k.  Ub  and  the  factor  3  represents  the  "  e-folding  ”  of  the  amplitude 
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after  the  onset  of  nonlinear  behavior.  Typically  this  predicts  values  for  (e8<j>)/Te  ^  1  which 
restricts  the  choice  for  the  electron  model.  In  particular,  a  linearized  electron  response  cannot 
always  be  expected  to  yield  correct  results. 

However,  if  the  beam  fraction  is  large  enough  and  the  beam  temperature  is  low  enough  so 
that  the  phase  velocity  lies  substantially  outside  the  ion  distribution  function,  then  the  per¬ 
turbed  saturation  potential  at  saturation  must  be  larger  than  predicted  by  Eq.  (13).  It  is  for 
these  cases  of  a  cold  and  strong  beam  that  ion  trapping  is  not  the  only  significant  saturation 
mechanism.  The  nonlinear  electron  response  and  displacements  of  the  background  ions  of 
x/Ln  ^1  play  major  roles  in  limiting  the  perturbed  potential.  Under  these  conditions  where  the 
background  displacement  becomes  nonlinear,  the  local  approximation  becomes  more  tenuous 
A  more  complete  nonlocal  treatment  will  be  considered  separately  in  a  later  publication. 

Closely  associated  with  particle  trapping  in  a  single  wave  is  an  enhanced  cross-field  motion 
of  the  particle  guiding  center  parallel  to  the  equilibrium  density  gradient.  As  shown  by  Nevins14 
this  guiding  center  constant  may  be  expressed  as 

vr  ~  *  constant  (14) 

where  the  field  has  the  form  <£  *=  <f>0(x)s\n(kyy  +  kzz  -  cut).  Therefore  large  changes  in  v.  are 
accompanied  by  large  excursions  parallel  to  the  equilibrium  density  gradient.  This  effect  is  most 
prominent  for  resonant  particles  since  resonant  particles  suffer  the  largest  changes  in  v- . 

The  electrostatic  field  is  polarized  almost  perpendicular  to  the  magnetic  field  and  so  it  is 
natural  to  consider  the  possibility  that  the  instability  may  be  saturated  by  perpendicular  trap¬ 
ping,  that  is  that  the  perturbed  electric  force  should  be  comparable  to  the  magnetic  force.  Set¬ 
ting  those  two  forces  equal  gives 


Comparing  Eq.  (13)  and  Eq.  (15) 
parallel  trapping  when 
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suggests  that  perpendicular  trapping  should  dominate  over 
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is  satisfied.  This  condition  is  very  difficult  to  satisfy  as  kj>s  =  1  and  y  <  0.lk,Ub  make  it  pos¬ 
sible  to  rewrite  Eq(15)  as 

Ub 

~  >  100  •  (17) 

LS 

B.  Motion  with  Several  Waves 

The  presence  of  several  waves  complicates  the  analysis.  First  of  all,  even  for  a  single 
value  of  ki  there  exist  two  unstable  waves.  One  wave  has  positive  phase  velocity  along  the 
magnetic  field  (  positive  k.  )  and  the  other  wave  has  negative  phase  velocity  (  negative  k.  ) 
along  the  magnetic  field.  The  wave  with  positive  phase  velocity  is  associated  with  instability  in 
the  beam  with  positive  velocity  and  similarly  the  wave  with  negative  phase  velocity  is  associated 
with  instability  in  the  beam  with  negative  velocity. 

However,  the  presence  of  two  waves  (  one  with  a  negative  phase  velocity  along  the  mag¬ 
netic  field  and  the  othe  with  a  positive  phase  velocity  along  the  magnetic  field  )  causes  no 
significant  change  over  the  case  of  a  single  wave.  The  reason  is  that  in  the  frame  of  one  of  the 
beams,  the  two  phase  velocities  of  the  unstable  waves  are  vastly  different.  The  wave  in  approx¬ 
imate  resonance  with  a  given  beam  is  responsible  for  nonadiabatic  change  of  that  beam.  The 
other  wave  is  of  very  high  frequency  and  only  causes  high  frequency  adiabatic  motion  of  the 
particles  in  that  beam. 

Next,  we  consider  the  case  of  multiple  values  of  k.  Frequently,  it  is  claimed  that  multiple 
waves  lead  to  quasilinear  diffusion.  However,  this  is  subject  to  two  constraints;  (1)  that  the 
wave  autocorrelation  time  is  short  compared  to  the  trapping  frequency  and  (2)  that  the  number 
of  unstable  waves  present  is  large.  Both  of  these  requirements  may  not  be  met  for  experimen¬ 
tal  or  for  the  simulations.  One  reason  for  this  is  that  all  of  the  unstable  waves  have  parallel 
phase  velocities  of  ~  ±  Ub  as  seen  in  Fig.  11.  Thus  the  wave  autocorrelation  time  is  very 
long.  In  addition  both  in  experiment  and  in  simulation  the  k  spectrum  is  discrete.  The  discrete 
nature  of  the  k  spectrum  is  an  important  consideration  in  most  simulations.  Typically,  in  a 
simulation,  a  single  mode  dominates  the  saturation  process. 
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A  more  complete  treatment  would  have  to  include  both  axial  and  radial  inhomogeneities. 
These  two  effects  may  alter  the  picture  given  in  the  preceding  paragraph  using  the  infinite 
medium  model.  The  complete  radial  problem  is  a  boundary  value  problem  and  as  such  may 
have  a  number  of  unstable  eigenfrequencies  for  a  given  value  of  ky  and  kz.  In  particular,  for 
those  cases  where  the  local  approximation  is  satisfied  the  eigenfrequencies  are  closely  spaced. 
Thus  the  complete  radial  problem  may  have  many  more  waves  with  comparable  growth  rates 
than  local  analysis  yields.  This  fact  should  tend  to  increase  quasilinear  effects. 

Axial  inhomogeneities  change  the  problem  by  causing  the  normal  modes  to  be  composed 
of  many  values  for  kz  rather  than  one  value  as  for  the  homogeneous  case.  Thus,  axial  variation 
of  equilibrium  parameters  may  assist  is  making  quasilinear  diffusion  a  more  important  mechan¬ 
ism  than  predicted  by  infinite  medium  theory. 

C.  Weak  Mode  Coupling 

Experimentally  it  was  seen  that  for  weak  cold  beams  there  were  peaks  in  the  spectrum  of 
the  noise  at  harmonics  of  the  ion  bounce  frequency  in  the  magnetic  well.  The  number  of  har¬ 
monics  seen  was  much  larger  than  the  likely  number  of  unstable  modes.  This  implies  a  possi¬ 
ble  mode  coupling  mechanism. 

The  theoretical  treatment  is  outlined  below  The  treatment  of  Sagdeev  and  Galeev4  is  fol¬ 
lowed.  The  particle  distribution  function  is  expanded  in  a  series  in  the  small  parameter  r\  as 

/« i/(iv  os) 

s  -  0 

and 

0=  £  0<*y  .  G9) 

J  -  0 

The  series  for  the  distribution  function  may  be  calculated  by  using 

S/(J)  -  T}  f  -  (20) 

where  the  integration  is  over  the  unperturbed  particle  trajectory.  Poisson’s  equation  then  closes 
the  system  of  equations  for  each  order  in  tj. 
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The  equation  giving  the  second  order  behavior  for  a  mode  (k,cu)  which  is  linearly  stable  is 
given  by 


-  -I 


e  (2) 


(k',k",6 


(i)  j  (i) 

<V  <V 


€  (k,<t>) 


(21) 


where  the  sum  is  over  all  modes  which  satisfy  the  conditions  k'  +  k  -  k  and  <o'  +  a/'  =  a>. 
This  expression  contains  a  very  complicated  function  e<2)  and  the  linear  response  given  by  e(1). 
Equation  (21)  actually  assumes  that  modes  (k,<o)  and  (k  ,to )  are  marginally  stable.  However, 
it  does  represent  the  correct  first  order  expansion  in  y/tu  and  so  retains  qualitative  justification 
even  for  finite  growth  rate.  One  relevant  fact  is  that  the  linear  response  may  become  small. 
This  is  realistic  if  both  the  beam  density  and  the  beam  thermal  spread  are  both  small.  Under 
these  conditions  harmonic  generation  is  possible. 

D.  Simulation  Results  in  the  Nonlinear  Regime 


As  indicated  above,  there  exist  several  nonlinear  mechanisms  associated  with  the  satura¬ 
tion  of  the  instability.  These  nonlinear  effects  will  be  demonstrated  be  using  the  results  for  the 
one  dimensional  code  and  the  two  dimensional  code.  In  addition  to  the  saturation  mechanisms, 
the  appearance  of  ion  phase  space  at  saturation  is  important  and  therefore  will  be  described. 

The  1-d  3-v  model  allows  only  one  of  the  beams  to  be  unstable  since  only  one  sign  of  ir¬ 
relative  to  kL  is  allowed.  Nevertheless,  this  model  is  a  useful  starting  point  for  examining  non¬ 
linear  effects.  For  example  Fig.  12  shows  particle  trapping  and  enhanced  cross  field  motion  near 
saturation  of  the  instability.  The  existence  of  the  guiding  center  constant  of  motion  dictates  the 
symmetry  seen  in  the  two  plots.  This  simulation  was  of  a  cold  strong  (  8  =  0.2  )  beam  with  the 
parameters  of  Sec.  3B.  and  as  a  result  nonlinear  electron  effects  and  nonlinear  motion  of  the 
background  ions  (  x/Ln  <  1  )  were  also  important  factors  in  saturation  of  the  instability. 

A  more  complete  model  of  the  instability  is  obtained  with  the  2-d  3-v  simulation  code, 
since  for  this  code  both  beams  are  unstable  in  the  simulation.  This  code  was  used  to  generate  a 
series  of  simulations  with  varying  thermal  spread.  The  thermal  spread  plays  a  large  role  in 
determining  the  magnitude  of  (eb<b/Te)  ,  at  saturation  and  the  appearance  of /(v.)  at  satura- 
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tion. 

The  simulation  parameters  are  the  same  as  used  in  Sec.  3B.  Table  1  gives  the  amplitude  of 
the  perturbed  potential  at  field  energy  saturation  Notice  that  the  Boltzmann  electron  response 
is  important  since  eh<t>/Te  is  not  neccesarrily  a  small  parameter.  Simulations  were  performed 
using  a  linearized  electron  response  and  the  saturation  levels  for  the  cold  beam  cases  were 
significantly  (  and  unphysically  )  enhanced. 

The  distribution  function  at  saturation  is  also  useful  for  demonstration  of  nonlinear 
effects.  These  distribution  functions  are  given  in  Fig.  13.  For  the  cooler  cases  the  entire  beams 
become  nonlinear  and  resonant  particle  trapping  is  not  seen.  As  the  beam  thermal  spread 
becomes  larger  and  larger,  obvious  trapping  effects  are  seen.  This  trapping  is  confined  to  the 
resonant  particles  and  the  bulk  of  the  distribution  function  remains  linear  as  shown  in  Fig.  14. 

Finally,  several  multimode  simulations  were  studied  to  examine  the  effects  of  retaining 
additional  stable  and  unstable  wavevectors.  Keeping  additional  unstable  wavevectors  caused 
only  slight  differences  in  the  single  mode  cases  mentioned  above.  This  is  due  to  the  discrete 
nature  of  the  simulation  k  space,  only  a  small  number  of  modes  have  comparable  growth  rates. 
For  example,  Fig.  15  shows  the  dominant  nonlinear  behavior  for  a  simulation  in  which  many 
modes  were  linearly  unstable.  As  can  be  seen  (  Fig.  15  )  the  constant  of  motion  given  in  Eq. 
(14)  is  still  very  useful. 

Keeping  additional  stable  wavevectors  at  multiples  of  the  unstable  wavevectors  does  lead 
to  mode  coupling.  This  effect  is  most  prominent  for  weak  cold  beams  as  mentioned  in  Section 
4B.  One  example  is  given  in  Fig.  16.  The  extent  of  coupling  obtained  in  this  matter  can  not  be 
compared  quantitatively  to  the  experiments,  since  the  experimental  results  depend  on  many 
more  parameters  than  included  here. 

The  nonlinear  effects  are  summarized  in  Table  2.  For  the  case  most  relevant  to  experi¬ 
mental  operation,  that  of  a  strong  thermal  beam,  the  saturation  is  by  beam  trapping  and  weakly 
nonlinear  electrons.  For  beams  with  very  large  thermal  spreads  parallel  to  the  magnetic  field, 
the  growth  rate  becomes  very  small  and  the  electrons  remain  linear. 
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5.  Axially  Inhomogeneous  Simulations 

This  section  is  a  preliminary  study  of  the  effects  of  axial  inhomogeneity  on  both  the  linear 
and  nonlinear  evolution  of  the  instability.  It  is  shown  that  considerable  shortening  of  the  axial 
length  is  necessary  for  significant  decrease  of  the  field  amplitude  at  saturation. 

In  order  to  relax  the  square  well  approximation  of  the  previous  section  an  axially  non- 
uniform  electron  response  was  used.  The  model  electron  response  is  given  by 

ne(z)  **  nQ<txv\e<t>/Tue<,z)\  (22) 

where  the  function  Tne(z)  represents  the  parallel  electron  temperature  and  is  a  function  chosen 
to  model  the  thermal  barrier  cell.  The  most  important  feature  that  the  model  must  retain  is  that 
the  effective  electron  temperature  increases  toward  the  middle  of  the  thermal  barrier  ceil. 

The  model  used  is  not  entirely  self-consistent  because  the  magnetic  field  is  assumed  to  be 
axially  uniform  for  pushing  the  ions.  In  reality,  the  reason  for  the  effective  electron  tempera¬ 
ture  to  vary  axially  is  that  energetic  electrons  are  trapped  in  an  axial  magnetic  well.  The  effect 
of  the  axially  varying  magnetic  field  can  be  neglected  for  the  ions  in  the  center  of  the  magnetic 
well  if  the  mirror  ratio  is  large  and  the  ions  are  injected  near  the  mirror  throat.  The  ion  motion 
in  the  stable  regions  then  just  gives  a  time  delay  and  the  details  of  that  motion  are  irrelevant. 

This  approach  makes  it  possible  to  examine  qualitatively  new  physics  as  compared  to  the 
square  well  model.  Unstable  regions  are  now  limited  to  the  region  in  the  center  of  the  cell 
where  the  instability  condition  U^k.  =  cs/(2L„ )  can  be  satisfied.  Further  away  from  the  well 
center  c5  is  much  smaller  and  unstable  values  of  k.  are  not  able  to  fit  into  the  system.  A  com¬ 
plete  theoretical  treatment  of  the  problem  must  include  all  of  the  bounce  resonances  in  the  well 
as  in  Sharp'5.  This  will  not  be  attempted  in  this  paper. 

In  order  to  perform  simulations  a  choice  must  be  made  for  the  function  Tn,(z).  Several 
different  functions  were  tried;  the  qualitative  effects  were  found  to  be  independent  of  the  exact 
functional  form  employed.  For  the  simulations  presented  here  the  functional  form  was  chosen 
(1)  to  be  symmetric  about  the  center  of  the  simulation  region,  (2)  to  be  constant  near  the 
edges  of  the  simulation  region,  (3)  to  be  parabolically  increasing  near  the  center  of  the 
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simulation  region,  and  (4)  to  be  a  continuous  function.  The  simulations  used  N.  —  65  with  the 
parabolic  region  confined  to  the  girds  larger  than  grid  number  8  and  smaller  than  grid  number 
57.  The  other  simulation  parameters  were  Ny  -  32,  6  -  0.2  cumax/fln  -  0.03,  v„/Ub  =  0.13 
and  Ub/cs  =  3.0  where  cs  is  evaluated  at  the  center  of  the  simulation  region.  Only  one  Fourier 
mode  was  kept  in  the  y  direction.  Approximately  25,000  particles  were  used.  The  system 
length  was  varied. 

Figure  17  shows  the  effect  of  axial  inhomogeneity  on  a  typical  potential  profile.  (  case  2  of 
Table  3  )  The  perturbed  potential  is  substantially  limited  to  the  center  of  the  simulation  region 
where  the  electron  Debye  shielding  is  the  weakest.  However,  the  potential  need  not  be  sym¬ 
metric  about  the  simulation  midplane  even  though  the  electron  response  is  symmetric.  Also, 
the  general  case  will  have  more  than  one  normal  mode  present  at  the  same  time  if  the  excita¬ 
tion  is  random  as  in  the  simulation  presented  here.  This  is  seen  to  be  the  case  from  Fig.  18 
where  a  time  history  of  the  real  part  of  the  potential  is  given  for  a  location  near  the  center  of 
the  system. 

The  saturation  of  the  instability  is  very  similar  to  that  of  the  square  well  simulations.  This 
can  be  seen  from  the  test  particle  trajectories  presented  in  Fig.  19.  Again  the  qualitative 
behavior  is  similar  to  that  dictated  by  the  adiabatic  constant  of  motion.  Table  3  shows  the 
values  of  the  perturbed  potential  at  saturation  as  the  system  length  is  varied.  Notice  that  a 
dramatic  decrease  is  not  attained  except  when  the  most  unstable  wavelength  is  longer  than 
about  twice  the  length  of  the  unstable  region. 

6.  Conclusions 

In  this  paper  a  simulation  model  has  been  used  to  examine  a  drift  wave  instability  which 
has  been  observed  in  some  experiments.  Qualitative  nonlinear  behavior  demonstrates  the 
importance  of  beam  trapping,  cross-field  diffusion,  and  nonlinear  electron  effects.  In  addition 
the  importance  of  the  thermal  spread  in  determining  the  saturation  behavior  has  been  shown. 
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Appendix  A. 


The  electrostatic  field  solver  used  in  the  simulations  is  described  in  this  appendix.  The 
field  solve  uses  either  (a)  linear  electron  Debye  shielding  in  conjunction  with  Poisson’s  equa¬ 
tion  (b)  quasineutral  Boltzmann  electrons  or  (c)  Boltzmann  electrons  in  conjunction  with 
Poisson’s  equation.  For  the  parameters  of  the  paper,  the  second  option  was  quite  adequate. 

The  field  solve  using  completely  nonlinear  Boltzmann  electrons  is  solved  for  in  a  manner 
using  techniques  given  in  Hockney  and  Eastwood16.  Given  an  initial  approximation  to  Poisson’s 
equation  df,  Poisson’s  equation  is  linearized  with  respect  to  that  approximate  solution  yielding 
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where  d>p~]  is  the  next  approximation  for  the  potential.  Subtracting  V2d>p  from  both  sides  of 


the  equation  yields  an  equation  for  the  small  quantity  W  =  <^+l— d>p  given  by 
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To  solve  this  equation  the  iterative  method  of  Concus  and  Golub17  is  used.  This  consists  of 
removing  (  approximately  )  the  W  depend  of  the  right  hand  side  of  Eq.  (A2).  This  yields 
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Here  AT  is  a  parameter  which  may  be  varied  to  improve  convergence.  The  optimum  value  is 
that  value  which  approximately  removes  the  IV  dependence  of  the  rhs.  Given  the  updated 
potential,  the  entire  procedure  is  started  again. 
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Saturation  Dependence  ot  Thermal  Spread 

Vr//  Ub 

- -l * 

0.01 - 

0.07 

- orn — 

- 

0725 

Tmax/  wmax 

0T3 

0.12 

Tj;09'5 — 

COJ73 

"  0.052 

\6o(f)/  Te)  sat 

U35 

0J3 

0729 - 

~XT22 - 

~ 0T5 - 

TABLE  1.  Variation  of  saturation  with  the  parallel  thermal  spread  of  the  ion  beams.  For 
larger  thermal  spreads  the  saturation  levels  should  decrease  very  quickly  since  for  those 
cases  trapping  is  the  sole  saturation  mechanism. 
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Saturation  Characteristics 

Vti/Ub 

0 

small 

large 

weak 

ion  beam  trapping;  enhanced 

trapping  of  resonant  ions; 

cross-field  motion  of  trapped 

enhanced  cross-field  motion 

ions;  coherent  mode  coupling 

of  trapped  ions 

strong 

nonlinear  electrons;  non¬ 

trapping  of  resonant  ions; 

linear  motion  parallel  to  the 

enhanced  cross-field  motion 

density  gradient;  ion  beam 

of  trapped  ions;  weakly  non¬ 

trapping;  enhanced  cross  field 

linear  electrons 

motion  of  trapped  beam  ions 

TABLE  2.  Nonlinear  properties  as  a  function  of  parameter  space. 
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Sam 

ration  Levels  as  a  Function  of  System  Length 

* 

(c  8(f)/  Te  )  saturation 

10.0 

0.38 

5.0 

037 - 

2.5  . 

0.15 

TABLE  3.  Saturation  levels  as  a  function  of  system  length.  The  parameter  is  the  sys¬ 
tem  length  divided  by  half  of  the  most  unstable  wavelength  evaluated  locally  at  the  center 
of  the  system.  The  parameter  Te  represents  the  electron  temperature  in  the  center  of 


the  simulation  region. 


2  k,  L 


FIG.  1.  Schematic  solutions  to  the  dispersion  relation  in  the  weak  beam  limit 


kj_^s 


FIG.  2.  Rea!  frequencies  for  a  numerical  solution  for  a  typical  set  of  parameters;  8  «  .1, 
Ub/cs  —  4.0,  vtl/Ub  -  0.05  ,  vJUb  —  0.0,  and  k:Ln  —  .1.  The  peak  in  the  growth  rate 
occurs  jjust  after  the  intersection  of  the  real  frequencies. 
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FIG.  3.  (a)  Normalized  maximized  growth  rate  maximized  over  the  k  vector  as 

a  function  of  A  =  L„/Lb  and  thfc  corresponding  values  of  (b)  the  normalized  real  fre¬ 
quency  a>r/comax  ,  (c)  kj> s  ,  and  (d)  k.Ln.  The  scale  length  of  the  background  plasma,  L„ 
,  is  held  fixed  while  the  scale  length  of  the  beam,  Lb  is  varied. 
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FIG.  3b. 
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FIG.  4.  (a)  Normalized  maximum  growth  rate  maximized  over  the  k  vector  as  a  function 
of  vj  Ub  and  the  corresponding  values  for  (b)  the  normalized  real  frequency  co,/co™x-  (c) 
k1ps  and  (d)  k:Ln. 
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FIG.  4c. 
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FIG.  5.  (a)  Normalized  maximized  growth  rate  as  a  function  of  v„/Ub  and  the 
corresponding  values  for  (b)  the  normalized  real  frequency,  (c)  kj>s,  and  (d)  k.L„. 
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FIG.  6.  (a)  Normalized  maximum  growth  rate  as  a  function  of  8  1.  8  and  the 

corresponding  values  of  (b)  the  normalized  frequency  ,(c)  kj>5  ,  and  W)  k.L„. 
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FIG.  6b. 
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FIG.  7.  (a)  Normalized  maximum  growth  rate  as  a  function  of  Ub/cs  and  the  correspond¬ 
ing  values  of  (b)  the  normalized  real  frequency,  (c)  pskL,  and  (d)  k.Ln. 
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FIG.  9.  Time  histories  of  (a)  the  square  of  the  mode  amplitude  and  (b)  the  square  of 
the  real  portion  of  the  potential  for  a  given  position  in  z.  For  this  simulation 
vn/Ub  =  0.19.  Random  excitation  was  used  to  initialize  the  simulation. 
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FIG.  10.  Theoretical  maximum  growth  rates  given  by  the  curve  are  compared  with  the 
simulation  results  given  by  the  crosses.  The  approximate  error  in  the  measurements  is 
given  by  the  width  of  the  crosses. 


FIG.  11.  Variation  of  (a)  the  real  frequency  and  (b)  the  growth  rate  in  the  k.L„  -  kj>s 
plane.  The  parameters  are  8  -  0  2  ,  [vi]0/^i>  “  ^  an^  v"/k*  “  ^0  The  con‘ 

tour  interval  in  (a)  is  0.883  wmax  with  the  maximum  of  0  9936  wmax  occurring  to  the  right, 
the  contour  interval  in  (b)  is  0.00552  a>;ax  with  a  maximum  of  0.0552  o>mix. 


0.8 


0.6 

KiPs 

0.4 

0.2 


FIG.  lib. 


2 


FIG.  12.  Particle  motion  (a)  showing  both  of  the  ion  beams  and  the  background  ions  in 
the  z  -  v.  plane  and  (b)  the  unstable  ion  beam  in  the  z  —  x/L„  plane.  Initially  all  parti¬ 
cles  are  at  x  -  0.  Since  the  background  ion  density  increases  with  increasing  x  a  positive 
value  for  x  indicates  that  the  particle  has  moved  in  the  direction  of  increasing  background 
ion  density.  This  convention  is  also  used  for  Fig.  16  and  Fig.  19. 


z 


FIG.  12b. 


FIG.  13.  Ion  distribution  functions /(v.)  at  saturation  for  the  five  cases  presented  in  Fig. 
10. 
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FIG.  15.  Test  particle  phase  space  for  the  multimode  simulation  (a)  in  the  z  —  v.  plane 
and  (b)  in  the  z  -  x/Ln  plane.  Initially,  the  test  particles  had  a  higher  speed  than  the 
magnitude  of  the  wave  phase  velocity  along  the  z  axis. 
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FIG  16.  Mode  energies  for  modes  one  through  three  where  mode  1  is  the  only  linearly 
unstable  mode.  The  other  modes  are  multiples  of  this  unstable  fundamental  mode  and 
have  finite  amplitude  due  to  mode-coupling.  The  simulation  parameters  are  v„/ U„  -  0.0, 
6  -  0.01,  Ub/cs  -  3.0,  n„Ar  -  0.40  ,  4096  particles  for  each  beam  and  the  background, 
and  there  are  65  grids  in  the  y  direction  and  33  grids  in  the  z  direction.  Only  three 
Fourier  components  of  the  potential  were  self  consistently  retained  in  each  direction. 


FIG.  17.  Magnitude  of  the  potential  squared  as  a  function  of  distance  along  the  magnetic 
field  for  a  time  where  the  instability  is  still  linear. 
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FIG.  18.  Square  of  the  real  portion  of  the  potential  near  the  center  of  the  simulation 
region  as  a  function  of  time. 
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FIG.  19.  Test  particle  trajectories  for  the  axially  inhomogeneous  simulation,  (a)  in  the 
z  —  v.  plane  and  (b)  in  the  z  —  x/Ln  plane.  The  test  particles  were  loaded  with  a  higher 
speed  than  the  wave  phase  velocity  along  the  z  axis. 
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